Examen final Modules 4 et 5
DUBii 2021
Template rédigé par les formateurs des modules 4 et 5 et modifié par Gaëlle Lelandais
14 April, 2021
Consignes
Complétez ce document en remplissant les chunks vides pour écrire le code qui vous a permis de répondre à la question. Les réponses attendant un résultat chiffré ou une explication devront être insérés entre le balises html code. Par exemple pour répondre à la question suivante :
La bioinfo c'est : <code>MERVEILLEUX</code>.N’hésitez pas à commenter votre code, enrichier le rapport en y insérant des résultats ou des graphiques/images pour expliquer votre démarche. N’oubliez pas les bonnes pratiques pour une recherche reproductible ! Nous souhaitons à minima que l’analyse soit reproductible sur le cluster de l’IFB.
Introduction
Vous allez travailler sur des données de reséquençage d’un génome bactérien : Bacillus subtilis. Les données sont issues de cet article :
Analyses
Organisation de votre espace de travail
# Je travaille dans mon répertoire de travail qui est associé au projet "dubii2021", sur le cluster de l'IFB.
cd /shared/projects/dubii2021/glelandais/modules-4-5-evaluation/.
# Création des différents répertoires qui seront utilisés pour organiser les fichiers.
mkdir raw-data
mkdir quality-controls
mkdir clean-data
mkdir mappingTéléchargement des données brutes
Récupérez les fichiers FASTQ issus du run SRR10390685 grâce à l’outil sra-tools [1]
# Travail dans le répertoire dédié au stockage des données initiales.
cd raw-data
# Comme je vais utiliser les outils en dehors d'un script (copier/coller),
# il me semble que cette commande est importante, afin de réserver des ressurces sur le cluster.
salloc --cpus-per-task=10 --mem=1G
# Chargement du module
module load sra-tools
# Lancement du téléchargement des fichiers FASTQ
srun fasterq-dump -p --split-files SRR10390685
# --> Deux fichiers FASTQ ont été obtenus.Combien de reads sont présents dans les fichiers R1 et R2 ?
# Chaque read est composé de 4 lignes dans un fichier FASTQ. La première débute
# par le symbole @. Je compte donc ces lignes dans les fichiers FASTQ et j'obtiens le
# nombre de reads. A noter que la même valeur est attendue pour les deux fichiers (données paired-end).
grep -c "^@" SRR10390685_1.fastq
grep -c "^@" SRR10390685_2.fastqLes fichiers FASTQ contiennent 7066055 reads.
Téléchargez le génome de référence de la souche ASM904v1 de Bacillus subtilis disponible à cette adresse
# Utilisation de la commande wget.
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/009/045/GCF_000009045.1_ASM904v1/GCF_000009045.1_ASM904v1_genomic.fna.gzQuelle est la taille de ce génome ?
# Pour simplifier les prochaines commandes, je décide de décompresser le fichier génome.
# A noter que étape n'est pas indispensable (ni même conseillée si l'on a le souci de préserver la quantité d'espace disque disponible).
gunzip GCF_000009045.1_ASM904v1_genomic.fna.gz
# Décompte du nombre de bases dans le fichier (attention de ne pas compter les caractères de fin de ligne !).
cat GCF_000009045.1_ASM904v1_genomic.fna | grep -v "NC_" | tr --delete "\n" | wc -c
# Note : J'ai trouvé de l'aide ici :
# https://dridk.me/genome_chiffre_1.htmlLa taille de ce génome est de 4215606 paires de bases.
Téléchargez l’annotation de la souche ASM904v1 de Bacillus subtilis disponible à cette adresse
# Utilisation de la commande wget.
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/009/045/GCF_000009045.1_ASM904v1/GCF_000009045.1_ASM904v1_genomic.gff.gz
# A nouveau, je fais le choix de décompression le fichier (étape non indispensable).
gunzip GCF_000009045.1_ASM904v1_genomic.gff.gzCombien de gènes sont connus pour ce génome ?
# J'écris une commande Bash qui permet d'extraire la troisième colonne du fichier et
# de compter des occurences du mot "gene".
cut -f 3 GCF_000009045.1_ASM904v1_genomic.gff | grep -c gene
# Remarque (après le rendu du projet) :
# --> Attention, après comparaison des résultats avec les autres stagiaires, je comprends que cette commande récupère aussi
# les "pseudogenes". Ainsi il aurait été préférable (pour avoir la bonne réponse), d'utiliser la commande awk incluant $3=="gene".4536 (incluant 88 pseudogenes) gènes sont recensés dans le fichier d’annotation.
Contrôle qualité
Lancez l’outil fastqc [2] dédié à l’analyse de la qualité des bases issues d’un séquençage haut-débit
# Changement de répertoire de travail.
cd ../quality-controls
# Chargement du module
module load fastqc
# Lancement des analyses qualités. Les fichiers de résultats sont écris dans le répertoire
# de travail.
srun fastqc ../raw-data/SRR10390685_1.fastq -o .
srun fastqc ../raw-data/SRR10390685_2.fastq -o .
# Création d'un rapport MultiQC
module load multiqc
srun multiqc -d . -o .La qualité des bases vous paraît-elle satisfaisante ? Pourquoi ?
- Oui
- Non
car les valeurs de score qualité sont élevées (> Q30) comme le montre le graphique “Per base sequence quality.” Les autres graphiques ne révèlent pas de problème particulier ici. Ce jeu de données est bien.
Lien vers le rapport MulitQC
Est-ce que les reads déposés ont subi une étape de nettoyage avant d’être déposés ? Pourquoi ?
- Oui
- Non
car tous les reads n’ont pas la même longueur. Un filtrage des séquences des adaptateurs a certainement été réalisé.
Quelle est la profondeur de séquençage (calculée par rapport à la taille du génome de référence) ?
# Je ne vais pas de solution pour répondre à cette question autrement qu'en réalisant
# le calcul moi même. Utilisation donc de la formule du cours :
# N*L/G ; avec N = Nombre de lectures, L = Longueur des lectures et G = Taille du génome
# ((7066055 * 2) * 151) / 4215606 = 506.2021 La profondeur de séquençage est de : 500 X.
Nettoyage des reads
Vous voulez maintenant nettoyer un peu vos lectures. Choisissez les paramètres de fastp [3] qui vous semblent adéquats et justifiez-les.
# Changement du répertoire de travail
cd ../clean-data
# Chargement du logiciel
module load fastp
# Allocation de ressources sur le cluster (j'ai augmenté ici la taille de la mépoire, car
# 1G n'était pas suffisant lors de mon premier essai).
salloc --cpus-per-task=10 --mem=5G
# Lancement de fastp. J'utilise par simplicité les mêmes paramètres que ceux utilisés lors du TP.
# La qualité des fichiers FASTQ est cohérente avec celle que nous avions pour les fichiers
# utilisés lors du TP.
srun fastp --in1 ../raw-data/SRR10390685_1.fastq --in2 ../raw-data/SRR10390685_2.fastq --out1 ./SRR10390685_1_clean.fastq --out2 ./SRR10390685_2_clean.fastq --html ./fastp.html --cut_mean_quality 30 --cut_window_size 8 --length_required 100 --cut_tail --json ./fastp.jsonLes paramètres suivants ont été choisis :
| Parametre | Valeur | Explication |
|---|---|---|
| cut_mean_quality | 30 | Valeur utilisée lors de la séance de TP. En dessus de cette valeur (30) la probabilité d’erreur pour une base est inférieure à 1/1000, ce seuil est classiquement utilisé. |
| cut_window_size | 8 | Valeur utilisée lors de la séance de TP. Elle indique le nombre de bases utilisées pour calculer une qualité moyenne. Cette valeur ne doit être ni trop faible, ni trop élevée. |
| length_required | 100 | Valeur utilisée lors de la séance de TP. Les lectures qui seraient plus petites que cette taille sont éliminées. |
| cut_tail | Utilisée lors de la séance de TP. La qualité des lectures est contrôlée égaklement du coté 3’. |
Ces paramètres ont permis de conserver 6777048 reads pairés, soit une perte de 4.1% des reads bruts.
Alignement des reads sur le génome de référence
Maintenant, vous allez aligner ces reads nettoyés sur le génome de référence à l’aide de bwa [4] et samtools [5].
# Changement de répertoire de travail
cd ../mapping/
# Chargement du programme
module load bwa
# Première étape de création des index
srun bwa index ../raw-data/GCF_000009045.1_ASM904v1_genomic.fna
# Deuxième étape des alignements des reads sur le génome
srun bwa mem ../raw-data/GCF_000009045.1_ASM904v1_genomic.fna ../clean-data/SRR10390685_1_clean.fastq ../clean-data/SRR10390685_2_clean.fastq > SRR10390685.sam
# Troisième étape de conversion des fichiers.
# Utilisation de samtools
module load samtools
# Succession d'étapes classiques (conversion en fichier binaire, tri et indexage).
srun samtools view SRR10390685.sam -b > SRR10390685.bam
srun samtools sort SRR10390685.bam -o SRR10390685.sort.bam
srun samtools index SRR10390685.sort.bam
# Pour terminer suppression des fichiers inutiles.
rm -f SRR10390685.sam SRR10390685.bamCombien de reads ne sont pas mappés ?
# Statistiques en relation avec les résultats de l'alignement.
srun samtools idxstats SRR10390685.sort.bam > SRR10390685.sort.bam.idxstats
srun samtools flagstat SRR10390685.sort.bam > SRR10390685.sort.bam.flagstat
# Affichage du contenu du fichier "flagstats"
less SRR10390685.sort.bam.flagstat
#13571369 + 0 in total (QC-passed reads + QC-failed reads)
#0 + 0 secondary
#17273 + 0 supplementary
#0 + 0 duplicates
#12826829 + 0 mapped (94.51% : N/A)
#13554096 + 0 paired in sequencing
#6777048 + 0 read1
#6777048 + 0 read2
#12746440 + 0 properly paired (94.04% : N/A)
#12769290 + 0 with itself and mate mapped
#40266 + 0 singletons (0.30% : N/A)
#0 + 0 with mate mapped to a different chr
#0 + 0 with mate mapped to a different chr (mapQ>=5)
# Le nombre de lectures non alignées est :
# 13571369 - 12826829 - 40266 = 704274
# Note : J'ai trouvé de l'aide ici :
# http://onetipperday.sterding.com/2020/11/explaination-of-samtools-flagstat-output.html704274 reads ne sont pas mappés.
Croisement de données
Calculez le nombre de reads qui chevauchent avec au moins 50% de leur longueur le gène trmNF grâce à l’outil bedtools [6]:
# Recherche des informations qui concernent le gène d'intérêt, création d'un nouveau fichier GFF
grep trmNF ../raw-data/GCF_000009045.1_ASM904v1_genomic.gff | awk '$3=="gene"' > trmNF_gene.gff3
# Chargement de la suite de programmes
module load bedtools
# Allocation de ressources sur le cluster
salloc --cpus-per-task=10 --mem=1G
# Sélection des reads alignés sur le gène d'intérêt.
# Il faudrait ajouter ici l'option -f pour indiquer que seuls les reads qui ont un chevauchement
# sur une longueur de 50% sont conservés.
srun bedtools intersect -a SRR10390685.sort.bam -b trmNF_gene.gff3 -sorted > SRR10390685_on_trmNF.bam
# Tri et indexage des reads
srun samtools sort SRR10390685_on_trmNF.bam -o SRR10390685_on_trmNF.sort.bam
srun samtools index SRR10390685_on_trmNF.sort.bam
# Nombre de reads présent dans le fichier résultat
srun samtools idxstats SRR10390685_on_trmNF.sort.bam > SRR10390685_on_trmNF.sort.bam.idxstats
# La solution précédente donne tous les reads qui chevauchent le gène d'intérêt. La condition 50% peut être ajoutée avec
# l'argument '-f 0.5'
srun bedtools intersect -a SRR10390685.sort.bam -b trmNF_gene.gff3 -sorted -f 0.5 > SRR10390685_on_trmNF_v2.bam
# Les autres traitements de données sont ensuite reproduits
srun samtools sort SRR10390685_on_trmNF_v2.bam -o SRR10390685_on_trmNF_v2.sort.bam
srun samtools index SRR10390685_on_trmNF_v2.sort.bam
# Nombre de reads présent dans le fichier résultat
srun samtools idxstats SRR10390685_on_trmNF_v2.sort.bam > SRR10390685_on_trmNF_v2.sort.bam.idxstats
less SRR10390685_on_trmNF_v2.sort.bam.idxstats2801 reads chevauchent le gène d’intérêt.
Visualisation
Utilisez IGV [7] sous sa version en ligne pour visualiser les alignements sur le gène. Faites une capture d’écran du gène entier.
Image finale
References
A work by Migale Bioinformatics Facility
https://migale.inrae.fr
Our two affiliations to cite us:
Université Paris-Saclay, INRAE, MaIAGE, 78350, Jouy-en-Josas, France
Université Paris-Saclay, INRAE, BioinfOmics, MIGALE bioinformatics facility, 78350, Jouy-en-Josas, France
